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A unique formulation of describing fluid motion is presented 
The method, referred to as "extended Lagrangian method," is inter- 
esting from both theoretical and numerical points of view. The 
formulation offers accuracy in numerical solution by avoiding nu- 
merical diffusion resulting from mixing of fluxes in the Eulerian 
description. The present method and the Arbitrary Lagrangian- 
Eulerian (ALE) method have a similarity in spirit— eliminating the 
cross-streamline numerical diffusion. For this purpose, we suggest 
a simple grid constraint condition and utilize an accurate discretiza- 
tion procedure. This grid constraint is only applied to the transverse 
cell face parallel to the local stream velocity, and hence our method 
for the steady state problems naturally reduces to the streamline- 
curvature method, without explicitly solving the steady streamline- 
coordinate equations formulated a priori. Unlike the Lagrangian 
method proposed by Loh and Hui which is valid only for steady 
supersonic flows, the present method is general and capable of 
treating subsonic flows and supersonic flows as well as unsteady 
flows, simply by invoking in the same code an appropriate grid 
constraint suggested in this paper. The approach is found to be 
robust and stable. It automatically adapts to flow features without 
resorting to clustering, thereby maintaining rather uniform grid 
spacing throughout and large time step. Moreover, the method is 
shown to resolve multi-dimensional discontinuities with a high level 
of accuracy, similar to that found in one-dimensional 

problems. © 1995 Academic Press, Inc. 


1. INTRODUCTION 

It is well known that fluid motion can be specified by either 
the Eulerian or Lagrangian description. Most CFD develop- 
ments over the last three decades have been based on the 
Eulerian description and considerable progress has been made. 
In particular, the upwind methods, inspired and guided by the 
work of Gudonov [1], have met with a great deal of success 
in solving fluid flows, especially where discontinuities exist. 
However, this shock capturing property has proven to be accu- 
rate only when the discontinuity is aligned with one of the grid 
lines since most upwind methods are strictly formulated in a 
one-dimensional framework and only formally extended to 
multi-dimensions. Consequently, the attractive property of crisp 
resolution of these discontinuities is lost. Even though research 
on genuine multi-dimensional approaches has recently been 
undertaken by several leading researchers, they are nevertheless 
still based on the Eulerian description. 


Recently, Loh and Hui [2] have convincingly demonstrated 
that a Lagrangian formulation can capture a contact discontinu- 
ity crisply, which is difficult to achieve by Eulerian formulation 
without resorting to some special treatment such as subcell 
resolution. Further developments have been carried out by Loh 
and Liou to solve real gas problems [3] and three-dimensional 
supersonic problems [4], The 3D extension is not so trivial, 
and in fact it involves somewhat tricky definition of cell move- 
ment and approximate solution of the multi-dimensional Rie- 
mann problem. Several interesting 3D problems that have not 
been attempted previously were calculated and again as in the 
2D case, high accuracy was achieved for resolving very com- 
plex shock-shock interactions. This method employs the point 
of view that it strictly follows the fluid particles released at 
some initial time line. The streamlines become a “time-like” 
coordinate and are used also for identifying particles. Therefore 
the method is naturally suitable for supersonic steady flow. No 
grid generation is needed a priori since the grid is part of the 
solution, viz., new grid lines are obtained as the solution marches 
in the time-like’ ’ direction. Unfortunately, restriction to super- 
sonic flows only limits the use of the method. To include the 
subsonic regime requires substantial conceptual changes. 

Numerically solving subsonic flows using the strict Lagran- 
gian concept becomes an excessive obstacle. Numerous re- 
searchers at Los Alamos, making substantial contributions to 
this subject, have proposed the so-called Arbitrary Lagrangian- 
Eulenan (ALE) method, for example [5, 6, 7] and others. Their 
method basically consists of two phases, namely Lagrangian 
and rezone phases. A third and implicit phase, similar to the 
Implicit Continuous-Fluid Eulerian (ICE) method [8] is inserted 
in the Lagrangian phase to allow for an efficient solution of 
flows at all speeds [5]. The method permits optimal use of grid 
which may move with the fluid (Lagrangian mode), or remain 
fixed with fluid moving through the cells (Eulerian mode), or 
move in any other prescribed manner (Arbitrary Lagrangian- 
Eulerian mode). This flexibility allows for the calculation of 
flows involving curved or moving boundaries and makes possi- 
ble calculations with minimum computational diffusion without 
excessive grid distortions. 

Although the concept for the Lagrangian approach and the 
equations to be solved are naturally very appealing, the proce- 
dure involving geometrical and variable constructs for the dis- 
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Crete system becomes rather complicated, see [5-7 1. Tracking 
the Lagrangian cells requires movement of cell vertices which 
must be done by interpolating variables and coordinates. It is 
well known that most multidimensional Lagrangian calcula- 
tions can only be continued for a finite time before the mesh 
is destroyed by “tangling,” or crossing of mesh lines [7]. The 
rezone phase is then added to prevent grid lines from crossing 
each other or grid tangling. Loss of accuracy, because of ac- 
counting for exchange of materials among the surrounding cells, 
can be inflicted through the continuous geometry and fiow p 
interpolations in response to the fluid particle motion [7|. In 
this regard, the paper can be thought of contributing a concept 
of moving the grid for this rezone phase in a specific way 
so that numerical diffusion at the cross stream faces remains 
vanished, as in the case of the strict Lagrangian approach. A 
more detailed comparison of the present approach and the ALE 
method is deferred to Section 3. 

Physically fluid particles at subsonic speed seem to adjust 
to motion and to surrounding (geometric or physical) constraints 
quickly and graciously, in particular sensing the upstream-prop- 
agating influences. Thus, a key to the success of a numerical 
Lagrangian procedure lies in how to properly and instantane- 
ously feed these upstream- propagating waves to the particles, 
while tracking them. It is indeed a very challenging research 
topic that motivates us to begin this exploratory investigation. 
This paper presents the salient features of the method, referred 
to as “extended Lagrangian method.” For flows at all speed 
regimes including purely subsonic and mixed flows, we demon- 
strate the advantages of the method over the Eulerian descrip- 
tion, with focuses on important features commonly seen in 
compressible flows, such as shocks, expansion waves, slip sur- 
faces, and interactions among them. We summarize in the fol- 
lowing some interesting features of the present extended La- 
grangian method. Specific advantages of the Lagrangian 
approach over the Eulerian one are given in Section 2 and 
additional advantages of the present extended Lagrangian ap- 
proach in Section 4. 

L The solution adapts to the flow variations (smooth or 
sudden), notably shocks and contacts, and as such it can be 
regarded as an automatic procedure for solution adaptation. 

2. Unlike the conventional adaptation techniques, there is 
no need for clustering grid points near the discontinuities. Very 
uniform grid can be maintained and in fact can also achieve 
orthogonality easily by construction. Thus discretization accu- 
racy does not deteriorate. Since streamlines do not cross, grid 
singularity or negative volume certainly will not occur. 

3. As will be seen later, the shock capturing quality in 2D 
is comparable to that found in the ID problem. This suggests 
that the present approach can be viewed as an alternative to 
the current genuine multi-dimensional approach. 

4. The contact discontinuity can be resolved crisply, since 
it is a streamline and as such no numerical diffusion is intro- 
duced due to fluid crossing the cell face. 


The rest of the paper is organized as follows. In Section 2 we 
compare differences of Eulerian and Lagrangian descriptions of 
fluid motion, with an emphasis on the numerical aspects. Some 
current methods relevant to the scope of the present paper con- 
cerning the Lagrangian and streamline-coordinates approaches 
are commented on in Section 3. Section 4 outlines the key ele- 
ments in formulating the present extended Lagrangian method 
for solving subsonic flows by retaining the advantageous fea- 
tures of the Lagrangian approach. A detailed formulation is 
then given in Section 5. Section 6 describes the grid motion 
of the present “extended Lagrangian” method. In Section 7 we 
briefly include the discretization method to complete the numer- 
ical procedure. Finally, in Section 8 we demonstrate the advan- 
tages of the proposed method by displaying solutions of flows 
at all speed regimes, containing various interesting features. 

It is noted that for the following sections, especially Sections 
2 and 5, while some of the statements may be well known, 
thus repetitive, to some readers, they are presented primarily 
with different concerns and motivations in mind, namely numer- 
ical viewpoint. In this way, a systematic and self-contained 
description can be made clear and we opt for elaboration rather 
than conciseness. 

2. EULERIAN VS LAGRANGIAN DESCRIPTION 

By definition, the Eulerian description oberves at fixed loca- 
tions the flow properties as fluid particles pass by. This has a 
close relation to the conventional computation approach in that 
each fixed grid point can be thought of as an observing station — 
corresponding to probes in measurement. With this approach, 
the meshes are generated mainly based on the geometry con- 
straint, with little regard given to the motion of fluid. Naturally, 
the grid lines will seldom coincide with fluid path lines. Even 
when grid lines are clustered near high-gradient regions using 
conventional adaptation techniques, they are not aligned with 
the particle path. A good example is the shock wave along 
which grid lines are densely distributed, and with which stream- 
lines make a nonzero angle, since the fluid will always pass 
through, not along, the shock wave. The angle is usually oblique 
in multidimensional flows. Consequently, the Eulerian approach 
has several severe numerical effects on the solution accuracy: 

1. Fluid particles are free to cross the grid line, thereby 
bringing (convecting) with them numerical mixing and diffu- 
sion across the cell interface. 

2. This numerical diffusion is only associated with the error 
resulting from approximating the convective terms. 

3. A contact/slip or shear layer is smeared ever increasingly 
w ith time and distance, unless some detection and special treat- 
ment techniques are employed. See [91 for example. 

Items 1 and 2 may appear to be two sides of a coin, indeed 
they are closely related on the discrete level. Item 1 can be a 
purely physical process if one were able to represent the cell- 
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face fluxes precisely as if in the continuous (differential) level. 
Thus, item 3 is the consequence of 1, in which the convection 
is represented approximately by a numerical procedure. It then 
becomes clear that the numerical procedure with minimum 
diffusion can be designed so that either the free crossing (ex- 
change) of mass is entirely restricted or the convective flux of 
the continuous system is precisely duplicated at the cell face. 
The latter proposition would seem impossible to do since the 
idea of discretization is to assume a function by an approximate 
valid in a finite domain, in which the approximate would tend 
to the function as the domain vanishes. Clearly the choice is 
the former. 

In spite of numerical diffusion resulting from approximation 
of convective terms, the Eulerian description does offer conve- 
nience and simplicity both conceptually and geometrically. The 
grid can be constructed regularly, simply based on geometry 
constraints and with little attention paid to the flow features. 
To enhance accuracy, grid adaptation is often applied to regions 
of high gradients. However, the concept of adapting to high 
gradients inevitably results in a skewed and distorted grid. This 
feature in turn will become more troublesome as two or more 
high-gradient areas intersect. 1 

It is quite safe to say that during the last three decades CFD 
algorithm researchers have primarily concentrated on devel- 
oping better (more accurate/robust/efficient) ways to deal with 
the convective terms , which exist only in the Eulerian formula- 
tion. Consequently, much success has been achieved, perhaps 
to the point of near perfection and little return could be gained. 
Unfortunately, inaccuracy due to numerical diffusion (mixing) 
in forming the interface numerical fluxes still exists and be- 
comes more exaggerated in multi-dimensional problems. On 
the other hand, since the convective terms do not explicitly 
appear in the Lagrangian formulation, the numerical mixing 
automatically disappears in the flux evaluation, rendering the 
Lagrangian approach attractive with respect to this viewpoint. 
However, other technical barriers have surfaced and discour- 
aged researchers from further pursuing development of methods 
for this approach. In what follows, we shall detail the concept 
of the Lagrangian approach from the viewpoint of numerical 
solution. The differences between the two descriptions will 
then become obvious. 

The Lagrangian description, by definition, states the motion 
and properties of given fluid particles as they travel to different 
locations. In particular, since the particle path in steady flow 
coincides with the streamline, no fluid particles will cross the 
streamline. In other words, while staying in contact, neigh- 

The unstructured grid approach will also face the same inaccuracy issue 
in approximating the convective terms, so long as the Eulerian description is 
adopted. To reduce the numerical diffusion error, adaptive remeshing techniques 
can be employed in both structured and unstructured grids by adding more 
and more points, i.c., finer and finer meshes, near the high-gradient regions. 
This viewpoint in fact deals with item 2 in the above comment — reduction of 
the the cross-streamline diffusion by shrinking the grid size. However, the 
remeshing comes at the expenses of computation and data structure. 


boring streams will not mix via convection, except in the molec- 
ular level where the physical molecular diffusion takes place. 
Therefore, the following numerical consequences can be 
realized: 

1 . No numerical diffusion is introduced across the cell inter- 
face since the computational cells follow the streamline, leading 
to a crisp resolution of contact/slip surfaces 

2. Fluid particles change motion (direction/speed) only as 
warranted, e.g., as shock or expansion waves are encountered. 
In other words, streamlines will bend, converge, or diverge 
only as situations demand. As a result, the shock capturing 
quality in 2D is comparable to that found in the ID problem, 
suggesting the present approach as an alternative to the genuine 
multi-dimensional approach. 

3. This description gives a realistic depiction of flow behav- 
ior; cells of same y-index form a streamline that is identifiable 
with flow visualization. 

The notion of using a Lagrangian approach to describe flow 
is not new. In fact, the very essence of following fixed particles 
also presents mathematical complexity to the approach, thereby 
limiting its scope of success. With the help of the new Lagran- 
gian formulations, numerical solution can be as easily obtained 
as for the Eulerian approach, only with the additional distinct 
advantages as stated above. In the following, we shall first 
review some current numerical procedures based on the Lagran- 
gian approach, commenting about their strengths and weak- 
nesses. Then we will focus on the applicability to the more 
difficult problem, namely the subsonic regime. 

3. REVIEW OF RELEVANT APPROACHES 

The present method, as will become clear in what follow, 
share several common concepts with the ALE method and the 
streamline-coordinate approaches. Yet, it also differs in either 
concept or procedure from the above two. Thus it is useful to 
give a brief discussion on these methods. 

3.1. Lag ra ng ian/A LE Meth ods 

The Arbitrary Lagrangian-Eulerian Technique (ALE) perhaps 
is the most well-knowm Lagrangian formulation in use at the 
present time. The technique, initially conceived and developed 
at Los Alamos during the 1970s, was implemented in several 
production codes such as CAVEAT and KIVA [10, 1 1 ], and the 
others. For the purpose of discussion, we will briefly summarize 
some of the key features in this technique; for a complete 
description of the procedure, the reader is referred to Refs. [5, 
6, 10, 11). The ALE method, while evolving throughout the 
years, maintains some basic constructs. It uses the staggered 
grid on which velocity components are stored at the cell verti- 
ces, thus directly giving displacements of vertices, and other 
variables at the cell center. The method basically consists of 
two phases of numerical procedures, namely Lagrangian and 
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rezone phases, using the time splitting concept. A third and im- 
plicit phase, similar to the Implicit Continuous-Fluid Eulerian 
(ICE) method [8) is added in the Lagrangian phase to allow for 
efficient solution of flows at all speeds [5]. The method permits 
optimal use of grid, which may move with the fluid (Lagrangian 
mode), or remain fixed with fluid moving through the cells (Eu- 
lerian mode), or move in any other prescribed manner (Arbitrary 
Lagrangian-Eulerian, ALE, mode). This flexibility allows for the 
calculation of flows involving curved or moving boundaries and 
makes possible calculations with minimum computational diffu- 
sion without excessive grid distortions. 

The Lagrangian approach and the equations to be solved are 
indeed very appealing conceptually. Tracking the Lagrangian 
cells requires movement of cell vertices which in turn involves 
interpolation of geometrical and variables data. Premature tan- 
gling failure can result from numerical errors such as discretiza- 
tion errors and nonuniform mesh, the latter being the subject 
of investigation in [7]. To prevent this grid difficulty, the rezone 
phase is added to regularize the cell distortion in a prescribed 
manner — optimal choice for controlling the grid distortion 
would be problem dependent. Loss of accuracy, because of 
accounting for exchange of materials among the surrounding 
cells in the ALE mode, can be inflicted through the continuous 
geometry and flow interpolations in response to the fluid particle 
motion. As a result, spurious error produced by this procedure 
can feed to the cycle of grid irregularity and tangling, see [7]. 

Recently a new Lagrangian Formulation was proposed by 
Hui and Van Roessel [12]. The inviscid conservation laws are 
transformed by using stream functions and Lagrangian time as 
independent variables. The stream functions serve to identify 
particles, while Lagrangian time represents time-like coordi- 
nate. Under this formulation, geometry conservation is enforced 
and each cell is literally a fluid particle. The grid movement, 
although often distorted, is obtained as part of the solution in 
the time-like (equivalently space) marching process from the 
derived geometrical conservation laws. The method has been 
shown to allow for extremely sharp resolution of contact discon- 
tinuities by Loh and Hui [2] in 2D and Loh and Liou [4] in 
3D problems. Multi-dimensional discontinuities are resolved 
with the same level of accuracy as their one-dimensional coun- 
terparts, indicating that the Lagrangian formulation inherently 
includes multi-dimensional flow characteristics. However, a 
severe limitation restricts the validity of the formulation [2-4] 
to only supersonic flows because the formulation is based on 
the use of the time-like coordinate. Thus, extention to subsonic 
flows based on the same framework appears impossible. 

Also, when a contact/slip surface exists initially, the fluid 
particles on either side of it may start separating from each 
other and eventually lose contact, disallowing interactions of 
fluxes (pressure only) from these fluid cells. Loh and Hui [2] 
applied a quick fix by subdividing the Lagrangian cells in a 
manner that the “time” increment (i.e. the sub-cell size) is 
adjusted to permit the cells remain in contact. Now, the strict 
“Lagrangian” concept begins to be relaxed since the sub-cells 


are only part of the initial cells, albeit belonging to the same 
identity. A formal transformation from the “Lagrangian time” 
(r-form) to the “Lagrangian” distance (A-form) was suggested 
by Hui and Zhao [13]. Conceptually, this draws close similarity 
to the steady streamline-coordinate method [14-18] in which 
one of the coordinates is a streamline. In fact, the discrete 
equations written in finite volume form from both formulations 
are identical for steady supersonic flows. Naturally, differences 
among these methods eventually appear as choices for repre- 
senting fluxes are made. Hence, it is relevant to discuss some 
aspects of the so-called streamline-coordinate method. 

3.2. Streamline-Coordinate Method 2 

Both the entropy and total enthalpy are conserved along 
streamlines if the flow is steady, inviscid, and adiabatic. Thus, 
there are advantages in exploiting these exact properties in the 
analysis by choosing streamlines as one of the coordinates. In 
the classical analysis, the other coordinate are defined to be the 
normals to the streamlines, this system is commonly known as 
Streamline Curvature [19], or Intrinsic, or Natural Coordinates 
[20]. Thus, in a two-dimensional irrotational steady flow, this 
coordinate system comprises of coodinates coinciding with the 
streamlines and equipotential lines. Being as natural as the 
intrinsic coordinate system may be, reports about solutions 
based on these coordinates, however, have been scarce in the 
literature. The streamline curvature method was popular in 
engine industry for calculating subsonic flows in turbomachin- 
ery. Nonconservative steady differential equations valid along 
and normal to the streamline are used and solved numerically 
by iterative relaxation procedures [14-16]. Making use of the 
advances in the Eulerian approaches, Giles and Drela [17, 18] 
solved the discrete approximate equations of the steady stream- 
line-coordinate equations in conservation form. Each 2D quad- 
rilateral cell is defined such that there is no mass across two 
of the four sides. Hence the mass and energy fluxes in each 
cell along the streamline coordinate are particularly simple. 
The only contribution to the steady-state equations from the 
streamline faces comes from the pressure contribution to the 
momentum equations. An auxiliary pressure relation must be 
assumed in order to close the resulting system. The grid posi- 
tions are not known a priori, and must be solved as part of 
the solution. Supersonic and transonic flows with shocks were 
solved using the artificial compressibility for maintaining stabil- 
ity [18]. A common feature in all these developments is the 
use of steady-state equations which are explicitly written on 
the basis of the streamline coordinates. This does not necessarily 
offer the best choice for resolving practical flows from the 
numerical point of view. The streamlines may be tightly bun- 
dled with large curvature, as in the vortical flow, or may run 
into singularities, such as separation and reattachment points. 

2 1 thank one of the referees who kindly brought my attention to the literature 
on the streamline curvature method. This section is the result of surveying 
some past works on the subject. 
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In these situations, we speculate that the grid (geometrical) and 
flow variables, as closely coupled, could result in numerical 
difficulty, eventually terminating solution. 

As stated in the last section, the so-called A-form of Hui and 
Zhao [13] for the steady supersonic flow indeed falls in the 
category of the streamline-coordinate method. While a rigorous 
algebraic manipulations has been carried out in the differential 
equation level, the discrete equations could be derived much 
directly from the integral consideration, as done in [17, 18] 
and later in Section 5. In their formulation, the A-coordinate 
coincides with the streamline and thus only pressure acts on 
the faces along A. The grid is obtained like the other variables 
in a space-marching fashion. From the forgoing discussion, it 
is clear that the methods classified in Sections (3.1) and (3.2) 
are identical in the case of steady supersonic flows and share 
the same objective for eliminating the numerical diffusion asso- 
ciated with fluid crossing cell boundaries. Having said this, we 
note that significant differences between approaches in Sections 
(3.1) and (3.2) exist for arriving at this objective, these differ- 
ences mainly bom out of numerical consideration. In this paper, 
we present an approach that can be seen to combine advantages 
of the above two: numerical flexibility of (3.1) for moving the 
grid and vanishing numerical diffusion (convective error) across 
the streams. Other distinct features will also be brought out in 
Section (5.3). 

In what follows we will first give the basic ideas for extension 
for subsonic flows in the next section and then describe detailed 
steps in Section 5. 


case. The net result is that we retain the essential beauty of the 
Lagrangian description that introduces no or minimal numerical 
diffusion across streamlines, but also at the same time we avoid 
its difficulty in grid tangling. On the other hand, we extend the 
streamline-coordinate method by not solving equations written 
in streamline coordinates, which are not known a priori, and 
by not using the auxiliary pressure relation [17, 18]. 

5. FORMULATION 

To facilitate the description, let us first define the notation for 
the relevant variables in the 3D Euler equations. The physical 
variables in a phase space of dimension 5 are denoted by a 
boldface uppercase letter or column vector whose elements are 
denoted by lowercase letters. 


r\ 


r\ 

pu 


f Pu ] 

pv 

u,.= 

pv 

pw J 


1 P w j 

l pe, ) 


\ ph.l 


where e, = e + 0.5 (u 2 + ir + w 2 ) and h, = e, + pip . The 
geometrical vectors in physical (Cartesian) space of dimension 
3 are denoted by an overhead arrow The fluid velocity is 

V = ui + vj + wk , (2a) 


4. EXTENSION TO SUBSONIC FLOWS 

A key element in the subsonic flow is the existence of the 
upstream-propagating wave. Thus, the existence of a body lo- 
cated downstream is transmitted to the oncoming fluid particles 
via this wave so that the particles can change motion accord- 
ingly. This immediately implies that we must abandon the time- 
like formulation cited above since it is only suited for pure 
initial value problems, such as supersonic flow where no influ- 
ence comes from downstream. Next, we must also abandon the 
idea of following a fixed particle, at least for the steady flows. 
Alternatively, we consider the steady streamlines as a set of 
lines that are occupied by particles released at the same location, 
different times and yet treated indistinguishably. The upstream- 
propagating influence is felt through the downstream particles 
on the same streamline in order to satisfy the governing conser- 
vation equations and boundary conditions in question. By de- 
scribing fluid motion along streamlines, we allow fluids to 
maintain their identity without tracking each specified particle. 
This definition is of course broader than and is a sufficient 
condition to the Lagrangian description, which follows motion 
of fluids of specific identity. Consequently, the present method 
is termed extended Lagrangian method. In other words, our 
approach (which boils down to the constraint given later in Eq. 
(10)) includes the strict Lagrangian description as a special 


and the normal vector of the boundary surface of a control 
volume 


S = s x i + s y j + s z k . 


(2b) 


The inviscid fluxes in 3D physical space are compactly writ- 
ten as 


P \ 

/°\ 

pu ' 


' pi 1 

pv 

V + 

PJ 

pw 


pi . 

ph ( 1 

\ 0 / 


where P - 


/°\ 

Pi 

Pj 

P * j 
\ 0 / 


(3) 


The first term in F is the flux of U ( convected by the fluid 
velocity V and the second term simply the pressure flux 
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FIG. 1. Definition of a control volume and material volume, (a) Control 
volume, denoted by “*'\ may be nonstationary, e.g., on rotating frame, (b) 
Material volume, V h = V 


For the following discussion, it is instructive to review some 
basic concepts used to describe fluids. Even though some state- 
ments may seem redundant to some readers, we attempt to 
arrive at essential conclusions primarily from the numerical 
point of view and for the sake of self-containedness. The deriva- 
tion given below is systematic and, to my knowledge, is not 
found elsewhere. 

To make the presentation self contained, a slight preliminary 
is useful. It is understood that the fluid has been considered 
to be a continuum. A convenient concept within continuum 
mechanics for describing a fluid motion is that of control vol- 
ume . In Fig. 1, let ft*(r) be a moving volume with bounding 
surface dft*(t); the local boundary velocity is V h . The volume 
is arbitrary and in general need not be identified with either 
physical boundary or specific motion of the fluid in ft*. Such 
a volume is called control volume. A special type of control 
volume is called material volume, denoted by ft(t), consisting 
of a collection of matter of fixed identity enclosed by a material 
surface, denoted by dft(r), of which every point moves with 
the local fluid velocity V . (See also Fig. 1.) If the volume ft(/) 
is shrunk to a point, the resulting material volume is called a 
fluid particle. Consequently, the fluid properties of the fluid 



FIG. 2. Definition of an Extended Lagrangian volume: the cell boundary 
dilf- is parallel to the fluid velocity V 


particle can be described mathematically in terms of space 
and time. 

Under the assumption of the continuum mechanics, let \ 
any continuous function, such as the density. The Reynolds 
transport theorem [21] gives the time rate of change of the 
“content” \ °f 


-j f xdv= f f xVh 9 dS, (4) 

dt J n*u) A J n*ui dt J A 


where the element surface dS of ft* is moving with the velocity 
V b . Note that V h may vary over the surface rift*. Again, a 
special case is when the theorem is applied to the material 
volume ft(/) with V h - V. 

The conservation laws (neglecting viscous fluxes for simplic- 
ity, without loss of generality in describing the approach) can 
be conveniently expressed over an arbitrary control volume 
ft*(f) in integral form: 


~L L [U,(V- v„) + P h ]-dS 

dt J J 


with P/ = P + 


= 0 , 

/ 0 \ 
1 0 ' 

0 

0 

\pV h I 


(5) 


From the above equations, three fundamentally different ap- 
proaches can result. 

5.1 Eulerian Description 

The Eulerian description assumes that the observer stays 
stationary with respect to the chosen frame of reference (e.g., 
inertial system). This requires: 

V h = 0 and ft* # ft*(f). (6) 

That is, the control volume is fixed in time. 

With the application of the Reynolds transport theorem Eq. 
(4), Eq. (5) is reverted to the familiar integral form: 

( f ™dv+ f JU,V + P]*</S =0, (7) 

J n f dt J 

where the superscript “E” is used to denote Eulerian frame of 
reference. In the discrete version of the above equation, each 
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FIG. 3. Generating an Extended Lagrangian computational cell so that the 
cell boundary (r,+ Uj+l — r iJ+l ) is parallel to V 

cell represents a control volume and is not moving in time, 
even though the flow may be unsteady. Note that the case in 
which an observer is fixed to a non-inertial frame, e.g., on 
rotating machines, also belongs to the Eulerian description. 

5.2 Lagrangian Description 

According to the strict definition, the Lagrangian description 
requires that the volume 0*(f) move with the instantaneous 



FIG. 4. Mach contours for M x = 0.4, 20% bump flow; (a) Lagrangian 
solution (134 X 70 grid), shaded area indicates numerical wake, (b) Eulerian 
solution (266 X 138 grid), and (c) pressure contours, noting normal behavior 
in the wake region. (Contour levels: „ = 0.0, Af ma , = 0.8, A M = 0.02.) 


fluid velocity and be identified as the material volume 11(f) (see 
also Fig. 1). That is, 

V h = V for all Vf > 0 and H*(f) = fl(f). (8) 
And the conservation laws are simplified to 

U dv+\ P„'ds=0. (9) 

cltiim jMm 

Clearly, the pressure remains as the only contribution to the 
flux on the surface. This renders extremely simple calculation 
of the time-rate of change of U dv — involving only the pressure 
acting on the bounding surfaces, only ifLL(t) is known precisely. 
However, the trajectory of the vertices is the part that often 
causes difficulties, resulting in large deformation or irregularity. 
(See [5-7].) Nevertheless, this is a very intriguing idea that 
avoids the nonlinearity in the equation of motion, thus reducing 
many difficulties associated with the convective terms that exist 
only in the Eulerian viewpoint. Such nice properties unfortu- 
nately have not been able to outweigh the drawbacks (see 
Introduction) and gain favorable reception over the Eulerian ap- 
proach. 

In the following, we propose a new approach that retains 
essential advantages of the above two approaches. 

5.3 Extended Lagrangian Description 

Close investigation of the surface integral reveals that the 
convective term can be eliminated also by requiring the con- 
straint: 

(V — V h )*S = 0. (10) 

This is the most general condition to vanish numerical diffusion 
resulting from convection. It amounts to either of the state- 
ments: (a) alignment of cell surfaces with the fluid relative 
velocity (V - V h ) — reducing to streamline coordinates for 
steady-state flow where V h — 0, or (b) V = V h , which is the 
strict Lagrangian description. Since (b) is included in Eq. (10) 
as a special case, our method is thus termed the “ Extended 
Lagrangian Method”, and Eq. (10) is termed constraint “EL”. 

It should be kept in mind that the condition, Eq. (10), is also 
valid at any instant , thus applicable to unsteady flow problems. 
Equation (5) becomes, as the constraint EL is imposed, 

T f kl Wv+ f „ P rdS 

dtJto'hi) 

+ L l uv-Vb)v t + i\ws = 0. 

The control volume now is denoted by superscript “EL” to 
indicate the present description. The surface dI2 EL is comprised 
of two types, d(2 EL - <H2 EL U df2 EL , where d!2 EL coincides 
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with the instantaneous particle paths and represents 
inflow and outflow faces. 

Clearly, the present “extented Lagrangian method” c 
bines the quality of the above two approaches: the sec 
surface integral includes both convective and pressure tet 
as in the case of Eulerian approach; the first surface inter 
on the other hand, merely has the effect of pressure, as in 
case of Lagrangian approach. That is, see Fig. 2, 
on dfl^ L <- both convected and pressure fluxes 
on dOp L only pressure flux 
For steady flow (hence V h = 0), there is no need for liter; 
following particles because no variations of motion with ti 
appear among the particles on the same streamline. Theref 


the question whether we strictly follow particles of same iden- 
tity or not is irrelevant in the formulation. Indeed, following 
the streamlines, rather than particles, is the essence of the 
present approach and this rescues us from facing the difficulty 
ot other Lagrangian approaches. Substituting V h = () m 
Eq. (1 1), we get 

"7 [ H XJdv ~h f P • dS 
dt J n fcl J *qel r ao 

c -4 (12a) 

+ J.< fU - y + P L^ = 0, in a*. 

together with the constraint, 
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V*S=0, ondilp 1 . 


(12b) M 


The semi-discrete form, including the time-dependent term re- 
tained for iteration purpose, can be cast as: 


- f n Vdv+ L P-' 

dt J u" mts g. 


.+ 2 

/Erlliy 


F, • S, = 0. (13) 



Examination of the above equations reveals some interesting 
insight. The inbalance of pressure in two neighboring cells wit 
common interface boundary (e.g., S,j in) cd “ ses 1 e 

change of flow direction (i.e., V in U„) of the fluids un 
consideration as well as change of their volumes CMj)- In other 
words, the deformation and dilatation of the fluid can be de- 
scribed. Indeed, the Lagrangian grid includes multi-dimensiona 
information and suggests how the fluid volume distorts in the 
flow. This point of view makes the description of fluid mo 
intuitively simple and clear. Moreover, it results in a major 
benefit in the numerical solution because this formulation 
avoids any arbitrary (numerical) mixing of fluids which in turn 
introduces numerical diffusion in the solution, notably across 
the contact discontinuity or shear layer. This diffusion error is 
common in the Eulerian formulation in which a nonstationary 
contact discontinuity is smeared without bound as time/space 
is marched. Furthermore, the advantage of the present approach 
is also clearly shown in its capability for crisply resolving 
multi-dimensional shocks. 

We note that the steady form of Eq. ( 1 3) is what the stream - 
line-coordinate method solves, see, e.g., 1 17, 1H|. Since t 
equation is valid only in the streamline coordinate, streamlines 
(thus grid lines) must be compatible with and be part of the 
solution. As a result, the number of unknowns is increased, 
leading to substantial increase in computation effort. Further, 
the effect of grids on the solution during the iterative 
process is strongly felt, robustness is of concern. Thus it is 
important to have a reasonably good initial position for the 

streamlines [ 18]. , 

In the present extended Lagrangian method, we are only 
interested in requiring that the limit form of Eq. (5) under the 
constraint, Eq. (10)— namely Eqs. (11) or (13)— be satisfied 
at convergence to unsteady or steady solutions. During interme- 
diate steps or any other times, it is not required that Eqs. ( 11 ) 
or ( 13) be solved. In other words, the equation to be solved for 
the present method is Eq. (5). The grid as well as conservation 
equations are not assumed to be in the streamline coordinate 
and the fluxes have the same form on every face of a cell. The 
grid instead evolves during the iteration according to Eq. (10) 
when it is invoked, resulting in a very simple and negligible 
effort to do. Since the grid is only a result of the solution an 
there is enough flexibility for adjustment hence goodness of 
initial grid or robustness is not an issue. 
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FIG. 6. Mach contours for M, = 0.675. 10% bump flow; (a) Lagrangian 
solution (134 x 70 grid), (b) Eulerian solution (266 X 138 grid). (Contour 
levels: M„„ = 0 0. M = 1.6, AM = 0.04.) 


We now conclude: the grid constraint EL expressed in Eq^ 

( 10) is all it takes to produce the extended Lagrangian method 
from an existing Eulerian method. In other words, ALE method 
or any other Eulerian codes can adopt Eq. (10) to arrive at the 
present method, leaving the existing code intact. 

In practice, the convective flux does not vanish exactly be- 
cause of the interpolation approximation used in following t e 
constraint EL. However, the error, namely the deviation from 
“zero” measured by the ratio of the convective and pressure 
fluxes on the transverse faces, max,, || U, V • AS||/||P • AS||, is 
negligible — on the order of 10 7 for theM = 0.4 case, discussed 

later in Section 8. 

We turn now to describe how Eq. ( 1 0) is fulfilled numerically. 

6. MOTION OF LAGRANGIAN GRIDS 

An important integral part of the present method is the grid 
motion that follows the constraint Eq. (10) or (12b) for steady 
flows. Two basic settings can be chosen for defining the motion 
of computational cells, namely the motion of cell centers or 
cell vertices. With the former approach the cell vertices will 
be defined by the position of neighboring centers, vice versa 
for the latter. Since the constraint, Eq. (12b), is imposed on 
the cell boundary, it is consistent to determine the vertex motion 
instead. This is easily done with the velocity field known from 
the solution. The constraint Eq. (12b) is equivalent to the kine- 
matic condition on a streamline. 


dx _ dy __ dz 
u v w 


(14) 


Again we remind the reader that a more general constraint is 
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FIG. 7. Distributions of variables on top and bottom walls for M* = 0.675, 10% bump flow; lines are Eulerian solution (266 X 138 grid) and symbols 
are Lagrangian solution (134 X 70 grid) with A and O denoting top and bottom walls, respectively. 


U — u h u — v b w — w h 

where V h can be chosen in a prescribed way that the cell remain 
regular, thereby eliminating grid tangling. 

As flow variables are defined at cell centers, the velocity 
components at cell boundary must be defined by some interpola- 
tion procedure from surrounding cells. In the present report, 
we outline the general notion of grid movement and give a 
specific strategy showing how the grid is moved to meet the 
constraint for the test cases included in the paper. Let us con- 
sider the two-dimensional cell (/, j), shown in Fig. 3. Three- 


dimensional cells can be treated similarly. Assuming the cell 
boundary is described by a line segment , — r /j4 , ). Since 
the segment is a part of a streamline, Eq. (14) gives 

.Vi+ 1.>4 1 = 1 4- nHx i+}J+ 1 - x iJ¥ i ) ( 1 6) 

where 

_ v 
m = ~ ■ 
u 

Here we list some possibilities for evaluating (m, v ): 
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prescribed number of flow variables iterations. As for steady 
flow calculations, it is not necessary to adjust the grid so fre- 
quently while the flow is still evolving. A more thorough inves- 
tigation about the optimal number of iterations per Lagrangian 
grid motion, its sensitivity to flow condition, and the effect of 
the velocity-averaging formulas will be useful, but is beyond 
the scope of the present paper. At this juncture, we note that 
for the results included in this paper, we used the mid-point 
average formula, Eq. (17a). 

To complete the numerical solution procedure, we define the 
numerical fluxes by employing the newly developed upwind 
scheme AUSM, which is described in full detail in [22, 23]. 
In what follows, we shall see that this new splitting has a very 
interesting bearing with the present extended Lagrangian 
method. 


FIG. 8. Mach contours for M* = 1.4, 4% bump flow; (a) Lagrangian 
solution (134 x 70 grid), (b) Eulerian solution (266 X 138 grid). (Contour 
levels: M mn = 0.4, M„ m = 2.2, = 0.045.) 


(a) Mid- point average 


V=l [2 V;,J + 2V,, t| + V I+1 , + V mj+I + V, ,, + V,-,J + ,1. 

(17a) 

(b) Upstream average, assuming u,j > 0, V/, j y 


V = zlV lJ +V,jn + V l .. lj +V i „h1. 


(17b) 


There are two unknowns, (x, y) /H j+i t in Eq. (16). Another 
condition is needed to complete the system. In this report, we 
prescribe the value of x-coordinate for each /th grid line, i.e., 
x — constant lines. This condition provides simplicity, but also 
yields accuracy as will be shown later. Furthermore, specifica- 
tion of x-coordinate allows one to put fine grids to resolve 
geometry details, e.g., near high curvature region. 

When the constraint, Eq. (12b), is satisfied for all j th grid line, 
the conservation laws are basically solved in a one-dimensional 
stream tube, because there is no flow across the 7-grid lines. 
As a result, high accuracy is expected with this virtually one- 
dimensional problem, which is demonstrated by the sharp repre- 
sentation of oblique shocks, as seen in Section 8 to be as 
accurate as their counterpart in one-dimension. The reason is 
that the formulation in itself already inherits multi dimensional 
information via the deformation of grid lines (i.e., streamlines). 

The cost of arriving at the above constraint is negligible even 
if it is done at each iteration, because calculation of Eqs. (16) 
and (17) is all it needs. The grid motion can be predicted in 
phase with the evolution of the flow variables, or else for a 


7. THE AUSM UPWIND METHOD [22, 23] 

To illustrate the concept, it is sufficient to consider only 
the one-dimensional system. As a first step, by recognizing 
convection and pressure as two physically distinct (but coupled) 
processes, we split the flux in the form of Eq. (3). In other 
words, these two terms deserve separate treatments. Mathemati- 
cally, we propose to separately deal with the genuinely nonlin- 
ear ({u - a, u -t- a) pair) and linearly degenerate ( u ) fields. 


F = u 



+ P = F ( . + P, 


F ( . = «U, 


(18) 


The overhead arrow “ has been dropped for we are con- 
cerned only with one-dimensional flow. Both Mach number 
and velocity splittings can be used to represent the convective 
quantity u in F t . In most cases, there is virtually no difference 
between calculated results of the two splittings. As found in 
[23], the velocity splitting is more robust in calculating unsteady 
shock tube problems by allowing a larger time step at start. 
Now, the numerical convective flux at the interface (denoted 
by subscript 1/2) straddling the left (L) and right (R) states, is 
effectively written as 


F c ]/2 — W 1/2U ( IVR» (19) 

where u m is the interface convective velocity. Let u m be writ- 
ten as: 


Wl/2 = «L+«R- (20) 

Several formulas are appropriate to define u~, e.g., 

f (« ± |w|)/2, if \u\ ^ a, 

u — \ 

[ ±(w ± aY 14a, otherwise, 


( 21 ) 
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FIG. 9. Distributions of variables on top and bottom walls for M * — 1.4, 4% bump flow; lines are Eulerian solution (266 X 138 grid) and symbols are 
Lagrangian solution (134 X 70 grid) with A and O denoting top and bottom walls, respectively. 


where a is the speed of sound. The convectible variable vector splitting. A differentiable pair of the ‘ + 1 and ‘ — ’ components 
U t is then upwinded solely based on the sign of u [l2 , viz. 


have been found to be effective. 




(U) L , 

(U) R , 


if u irl s 0, 
otherwise. 


(22) 


P = 


| Ml ± sgn(M))/2, 

I p(M± \)\2 + M)/4, 


if |w| ^ a , 
otherwise. 


(24) 


We turn now to the pressure term by writing: 

Pm. = Pi + Pv. • (23) 

Similarly, a whole host of choices are possible for the pressure 


This completes the definition of the numerical flux F. Putting 
( 1 9) and (23) together, we recast the interface flux in the follow- 
ing form 


F 1/2 — W 1/2 2 (U,L + U ( -r) 2 l Wl/2 l 


(25) 
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FIG. 10. Mach contours for M* = 1.8, 15% ramp flow; (a) Lagrangian 
solution (100 X 41 grid), (b) Eulerian solution (198 X 80 grid). (Contour 
levels: M mn = 0.4, M™ = 2.2, AM = 0.045.) 


where A l/2 { • } = {*} R — { * } L ♦ Here the first term on the RHS 
is clearly not a simple average of the ‘L 1 and k R’ fluxes, but 
rather a weighted average via the convective velocity. The 
dissipation term has merely a scalar coefficient |w,^| and requires 
only O(n) operations for an ^-dimensional vector F. Further- 
more, since there is no differentiation (or Jacobian matrix) 
involved in evaluating F !/2 , the present method is easily ex- 
tended to a general equation of state and non-equilibrium flows 
and the cost is only linearly increased with the additional conser- 
vation equations considered. Unlike the Roe or Osher schemes, 
the extension does not yield additional ambiguity such as the 
definition of averaged or intermediate states. Also, numerical 
tests strongly suggest an entropy-satisfying property by the 
present method. 

To achieve higher-order spatial accuracy, a MUSCL-type 
procedure is followed to upwind-extrapolate variables (primi- 
tive variables in the calculations presented in this paper), with 
TVD limiters incorporated [24]. Then, a two-stage Runge- 
Kutta procedure is used to integrate the semi-discrete system 
Eq. (13), subject to the kinematic condition Eq. (12b). 

The subsonic inflow conditions are imposed by specifying 
total enthalpy, total pressure, and flow angle, while the outflow 
conditions are obtained with specified static pressure and ex- 
trapolated total enthalpy, total pressure, and flow angle. The 
usual tangency procedure is used at the cell boundary that 
coincides with a physical wall — no ghost cells are used. The 
wall pressure is gotten using linear extrapolation from interior 
data, and so is the total enthalpy. 

8. TEST PROBLEMS AND DISCUSSION 

We will show examples for flows at all speed regimes, featur- 
ing solution accuracy and grid aspects. The plots are organized 
uniformly for all cases presented. We show Mach contours 


overlaid on the grid used. Fine grids, by doubling grid number 
in both directions, are used only in the Eulerian calculation for 
comparison purpose. There are 80 and 160 cells on the bump 
respectively for the coarse and fine grids. The grids shown are 
only one half of the whole grid in the Lagrangian case, and 
one quarter in the fine-grid case, thus corresponding to roughly 
the same location in both plots. The symbols denote the en- 
tended Lagrangian solutions and the lines are the Eulerian 
solutions. The Mach contours are chosen for presentation so 
that any numerical anomaly or inaccuracy can be more easily 
depicted than the other variables such as pressure. 

The first example is the purely subsonic flow in which an 
M* = 0.4 flow enters a channal with 20% circular bump, see 
Fig. 5 for detailed geometry. The Mach contours, given in Fig. 
4, show nearly perfect symmetry about the midchord, except 
in the wake region. The wake region suggests an entropy pro- 
duction (numerical diffusion), likely due to numerical wall 
boundary condition, which to my knowledge is still a gray area 
in CFD. The Lagrangian solution appears to result in a narrower 
numerical wake, as indicated in the shaded area; the Lagrangian 
result is roughly half the width of the Eulerian result. Moreover, 
close examination of these plots reveals that the numerical 
wake of the Lagrangian solution does not spread while on the 
contrary the wake of the Eulerian solution grows with distance. 
It is worth noting that the present grid automatically evolves 
from initial Eulerian grid into the grid system seen in Fig. 4, 
according to Eqs. (16) and (17a). The grid spacing between 
streamlines increases near the stagnation points and converges 
as flow accelerates. The detailed distribution of variables on 
the top and bottom walls are plotted in Fig. 5. The fine-grid 
Eulerian solutions are included for comparison. The agreement 
is remarkable and again symmetry is quite evident. However, 
the fine-grid Eulerian solution over-predicts the pressure at the 
leading and trailing stagnation points, whose theoretical value 
is 1.116, while the Lagrangian solution closely matches — 
mainly due to the expansion and contraction of grid lines near 
these two points with the largest change in slope. Hereafter, 
for briefness and contrast to the Eulerian solution, we shall 
take the liberty of loosely using the term “Lagrangian solution 1 ’ 
to mean the solution obtained by the extended Lagrangian 
method described in this paper — not in the strict Lagrangian 
sense. It is also noted that the fine-grid solution took consider- 
ably more iterations than the Lagrangian solution to converge. 

The second example involves the popular test of transonic 
flow in a channel with 10% bump and A/* = 0.675. Results 
are given in Figs. 6 and 7. Again, we show the Mach contours 
and distributions on both walls. The agreement of the coarse - 
grid Lagrangian solutions with the fine-grid Eulerian solutions 
is excellent. The shock resolution from the Lagrangian method 
is outstanding, so is the prediction of the so-called “Zierep” 
singularity at the foot of the shock on the curved surface. The 
Lagrangian solution again yields a grid that senses the global 
flow characteristic. Notice that no clustering of grid is necessary 
to capture the shock in the correct location. 
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The next case is also a standard one, involving an = 1.4 
flow and 4% bump. This case consists of small subsonic pocket 
resulting from short Mach stems on both walls, shown in Fig. 
8. Shock-shock-expansion waves interactions take place be- 
hind the trailing edge. The shock locations are in excellent 
agreement with the fine-grid solutions. Close examination of 
the Lagrangian grid shows that the grid lines (also streamlines) 
remain straight until the shock is encountered, as it should. The 
Eulerian grid however already began bending at the leading 
edge for all j th lines, simply because of geometry constraint. 
The present method yields grids that are conforming with the 
flow features by bending, expanding/contracting. The net result 
is that excellent shock resolution is obtained, even though the 
shock is oblique to the grid line. Figure 9 displays a one-cell 
capturing of the oblique shock. This is not entirely surprising 
since the present formulation has already taken account of 
the multi-dimensional nature of the flow via the streamline 
deformation caused by the fluxes (pressure forces) of sur- 
rounding fluids. 

The fourth test is an = 1.8 flow over a 15° ramp. This 
case consists of a Mach stem about 10% of channel height, a 
slip line emanating from the triple point, and reflected shocks. 
In Fig. 10, the mach contours depict an overall picture of the 
flow, demonstrating a sharp resolution of the ramp shock, Mach 
stem, and the subsequent shocks. The slip line, whose strength 
is being weakened by the expansion wave generated at the 
ramp shoulder and transmitted through the first reflected shock, 
is resolved to the same level of accuracy as given by the fine- 
grid solution, i.e., with the same level of spatial spreading. 
Figure 1 1 vividly displays how the grid lines bend as the shock 
is encountered and change direction according to the flow. The 
grid itself already suggests the flow structures, train of shock 
reflections, expansions, as well as the Mach stem across which 
there is no change of flow angle. It is worth noting the clear 
slipline emanating from the triple point. In contrast to the shock- 
aligned grid, the present grid is aligned with the streamlines, 
which will never cross each other, but the shocks can. Thus 
the present method is indifferent to whether the high-gradient 
regions intersect. The profiles (Fig. 12) on the walls show 
good agreement of both solutions. Excellent shock resolution 
capability is observed on both walls even the second reflected 
shocks remain well resolved. 

CONCLUDING REMARKS 

We have presented a unique formulation for dealing with 
subsonic as well as supersonic flows. The method, referred to 
as extended Lagrangian method, combines the accuracy belong- 
ing to a Lagrangian description and the robustness and simplic- 
ity of an Eulerian description. Through systematic comparison 
with fine-grid solutions and a theoretical check for flows at 
various regimes, we have demonstrated its capability for crisply 
capturing high-gradient regions that are not aligned with the 
grid. In contrast to adaptive approaches reported to date, the 


present approach already automatically inherits the ability to 
adapt to flow features, but based on an entirely different adap- 
tive philosophy in which there is no need for clustering grid 
lines. Since the grid lines are not aligned with high-gradient 
areas, they are maintained regular and uniform. Moreover, one 
set of grid lines that coincides with streamlines depicts vividly 
a form of “numerical flow visualization.” Without resorting 
to arbitrary detecting criteria, the present approach not only 
predicts well the nonlinear waves, such as shock and rarefaction 
waves, but also is especially amenable to treating a linearly 
degenerate field, such as a contact discontinuity. Furthermore, 
since the grid spacing is maintained relatively unform, a large 
time step is permitted throughout the calculation, thus increas- 
ing efficiency. Also the common adaptive strategy will have 
great difficulty in the case of intersecting shocks, because the 
grid lines will cross each other, if not checked. We also suggest 
that the present extended Lagrangian method is a viable alterna- 
tive approach to the current multi-dimensional scheme and grid- 
enriching adaptive procedure for complex flows having high- 
gradient regions. Further development and 3D applications are 
currently underway and will be reported in the future. 
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